Write a custom CUDA kernel to optimize `torch.nn.GaussianNLLLoss`.

The original operation is defined by the formula:
`loss = 0.5 * (log(max(var, eps)) + (input - target)^2 / max(var, eps))`

**Problem Analysis:**
The standard PyTorch implementation is a textbook example of a memory-bound operation. It executes a long chain of element-wise CUDA kernels:
1. `clamp(var, min=eps)`
2. `subtract(input, target)`
3. `square()`
4. `log()`
5. `divide()`
6. `add()`
7. `multiply(0.5)`
Each of these steps materializes a full-sized intermediate tensor in global GPU memory, leading to extremely high, unnecessary memory bandwidth consumption. This is the primary performance bottleneck, followed by the overhead of launching multiple kernels and a final reduction pass.

**Optimization Strategy: Fully Fused Computation and Parallel Reduction**

The strategy is to fuse this entire chain of operations, including the final reduction, into a single, efficient CUDA kernel pass.

1.  **Complete Operator Fusion**: All seven mathematical operations are fused into a single calculation performed within each CUDA thread. The thread reads `input`, `target`, and `var` once from global memory, computes the final loss value using only registers, and writes the result. This completely eliminates all intermediate tensors. The `max(var, eps)` clamping for numerical stability is handled inside this fused operation.

2.  **Element-wise Fusion (`reduction='none'`)**: A single kernel is launched where each thread computes the full, fused loss formula for one element and writes the result directly to the output tensor.

3.  **Fused Parallel Reduction (`reduction='mean'/'sum'`)**: For reduction modes, the fused computation is combined with an efficient parallel reduction algorithm:
    *   **Stage 1 (Calculation & Block-Level Reduction)**: Each thread block is assigned a subset of the data. Threads compute the loss for their elements, and then collaboratively sum these losses using a fast, tree-based reduction in **shared memory**. Each block writes its single partial sum to a small temporary global buffer.
    *   **Stage 2 (Final Reduction)**: A simple `tensor.sum()` is called on the small buffer of partial sums to get the final total loss. For 'mean', this sum is then divided by the total number of elements.

This approach transforms the multi-stage, memory-intensive operation into a single-pass, compute-efficient kernel that maximizes GPU utilization and dramatically reduces execution time.
  
Here's an example to show you the syntax of inline embedding custom CUDA operators in torch: The example given architecture is:   
  
```python
import torch
import torch.nn as nn

BATCH_SIZE = 512
DIM = 4096
SHAPE = (BATCH_SIZE, DIM)
REDUCTION = 'mean'
EPS = 1e-6

class Model(nn.Module):
    def __init__(self, eps=1e-6, reduction='mean'):
        super(Model, self).__init__()
        self.loss_fn = nn.GaussianNLLLoss(full=False, eps=eps, reduction=reduction)
    
    def forward(self, input_tensor: torch.Tensor, target_tensor: torch.Tensor, var_tensor: torch.Tensor) -> torch.Tensor:
        return self.loss_fn(input_tensor, target_tensor, var_tensor)

def get_inputs():
    input_tensor = torch.randn(SHAPE, dtype=torch.float32)
    target_tensor = torch.randn(SHAPE, dtype=torch.float32)
    # var 必须是正数
    var_tensor = torch.rand(SHAPE, dtype=torch.float32) + EPS
    
    return [
        input_tensor.contiguous(), 
        target_tensor.contiguous(), 
        var_tensor.contiguous()
    ]

def get_init_inputs():
    return [EPS, REDUCTION]